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GAS-KINETIC THEORY BASED FLUX SPLITTING METHOD FOR IDEAL 

MAGNETOHYDRODYNAMICS 

KUN XU* 


Abstract. A gas- kinetic solver is developed for the ideal magnetohydrodynamics (MHD) equations. The 
new scheme is based on the direct splitting of the flux function of the MHD equations with the inclusion 
of “particle” collisions in the transport process. Consequently, the artificial dissipation in the new scheme 
is much reduced in comparison with the MHD Flux Vector Splitting Scheme. At the same time, the new 
scheme is compared with the well- developed Roc-type MHD solver. It is concluded that the kinetic MHD 
scheme is more robust and efficient than the Roe-type method, and the accuracy is competitive. In this paper 
the general principle of splitting the macroscopic flux function based on the gas- kinetic theory is presented. 
The flux construction strategy may shed some light on the possible modification of AUSM- and CUSP- type 
schemes for the compressible Euler equations, as well as to the development of new schemes for a non-strictly 
hyperbolic system. 
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1. Introduction. The development of numerical methods for the MHD equations has attracted much 
attention in the past years. Godunov- type schemes arc considered particularly useful here. On the basis 
of Roe’s method [18], Brio and Wu developed the first Flux Difference Splitting (FDS) scheme for MHD 
equations [3]. Aslan also followed the idea of fluctuation approach to construct a second-order upwind MHD 
solver [1]. Zachary et al. applied an operator splitting technique and devised a high-order Godunov type 
method [28]. During the same period, the multidimensional extension of MHD solvers was done by Ryu et al. 
[19] and Tanaka [22]. On the basis of the nonlinear Riemann solver, Dai and Woodward extended the PPM 
method [5]. Powell et al. constructed an eight- wave family eigensystem for the approximate Riemann solver 
[15, 16]. Most recently, based on the Lax-Friedrich flux splitting technique, Jiang and Wu applied a high-order 
WENO interpolation scheme to the MHD equations [9]. In order to increase the robustness and simplify the 
complicated Roe-type MHD solver, based on the HLL method, Linde developed an adequate Riemann solver 
for the heliosphere applications [10]. A majority of the methods mentioned above applied characteristic 
decomposition for the MHD waves, where the entropy, slow, Alfven and fast waves have to be considered in 
the evaluation of a single flux function. Because of the wave decomposition procedure, considerable work is 
required to evaluate and justify the MHD eigensystem, where the non-strictly hyperbolicity causes additional 
difficulty. Due to the same reason, the issue related to the direct extension of the Flux Vector Splitting (FVS) 
scheme to the MHD equations was hardly addressed. The search for robust, accurate and efficient MHD 
flow solvers is still one of the primary directions in MHD research. 

For the Euler and Navicr-Stokes equations, the development of gas- kinetic schemes has also attracted 
attention [25]. A particular strength of kinetic schemes lies precisely where Godunov- type FDS schemes 
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often fail, such as carbuncle phenomena, positivity, entropy condition [26]. However, like any other FVS 
method, the Kinetic Flux Vector Splitting (KFVS) scheme is very diffusive and less accurate in comparison 
with the Roe- type Riemann solver, especially for shear and contact waves. The diffusivity of the FVS 
schemes, such as Steger- Warming, van Leer and the KFVS [21, 24, 17], is mainly due to the particle or wave- 
free transport mechanism, which sets the CFL time step equal to particle collision time. Consequently, the 
artificial viscosity coefficient is always proportional to the time step. Even though numerically the high-order 
FVS methods can get crisp shock resolution by using a MUSCL-type reconstruction method, physically it 
is impossible to develop a second-order FVS scheme without correcting the free transport mechanism. In 
order to reduce the diffusivity, particle collisions have to be modeled and implemented in the gas evolution 
stage, such as that in the BGK scheme [27]. 

The construction of the gas-kinetic FVS scheme for the MHD equations started from Croisille et al. 
[4], where a MHD KFVS solver was obtained by simply extending the KFVS flux function of the Euler 
equations. The above MHD KFVS scheme is very robust and reliable, but over- diffusive, especially in the 
contact discontinuity regions [10]. Recently, another interesting gas-kinetic MHD solver has been developed 
by Huba and Lyon [6]. Different from the earlier approach, Huba and Lyon constructed two equilibrium 
states and a transport equation to recover the MHD equations. An important aspect of this method is that 
it provides a framework to incorporate additional terms into the MHD equations, e.g. anisotropic ion stress 
tensor, and anisotropic temperature distribution. However, the physical basis of the transport equation and 
the reliability of the equilibrium states need to be further investigated. Since Huba and Lyon’s flux function 
keeps the FVS nature, large numerical dissipation is intrinsically rooted. 

In this paper, we are going to construct a new kinetic flux splitting method for MHD equations. Based 
on the BGK-typc scheme, the KFVS MHD solver of Croisille et al. is generalized by including particle 
collisions. As a result, the new scheme reduces the numerica: dissipation significantly and gives a more 
accurate representation of wave interactions. In Section 3, the lew scheme is compared well with the Roe- 
type MHD solver [3, 16]. The flux construction method presented in this paper splits the macroscopic flux 
function directly, therefore it is very useful in the design of numerical methods for complicated hyperbolic 
system. 

2. Gas-Kinetic Approach to MHD Equations. In the one-dimensional case, the MHD equation 

qt + F{q) x = 0 

has the following forms [3] 

Pt + (pu)x — 0, 

(pU) t + (pU 2 + p, - B 2 ) x = 0, 

(pV) t + ( pUV - B x B y ) x = 0, 

(2-1) < ( P W) t + ( pUW - B X B Z ) X = 0, 

(B y ) t + (B y U - B X V) X = 0, 

(B z ) t + (B Z U - B X W) X = 0, 

, (pe) t + {(pe+p.)U - B X {B X U + B y V + B Z W) = 0, 
where p* is the total pressure 

p.=P+\(B 2 x + B 2 y + 3 2 z ), 
and p is the gas pressure. The total energy density is 

pe = \p{U 2 + V^ 2 + W 2 ) + pe+ 1(.3 \ + B 2 y + B 2 Z ). 
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For ideal equilibrium gas, the internal energy is related to pressure through the relation 

pe = p/{ 7 - !)• 

Due to having a different physical origin, it should be emphasized that in order to properly split the energy 
flux function, the splitting of internal energy flux peU and the splitting of work done by the pressure pU 
should be different, although they are only different by a constant 1/(7 — 1) for the ideal gases. 

Theoretically, it is very difficult to construct an equilibrium state and a transport equation to exactly 
recover the above ideal MHD equations. However, instead of constructing the equilibrium distribution for 
the flow and magnetic field, we can split the MHD flux function directly on the macroscopic level with the 
consideration of gas- kinetic theory. 

2.1. Gas-Kinetic Flux Splitting Method. In the gas-kinetic theory, the flux is associated with 
the particle motion across cell interface. For ID flow, such as the ^-direction, the particle motion in this 
dimension is most important in the determination of the flux function. Other quantities, such as y-direction 
velocity, thermal energy, and magnetic field, can be considered as passive scalars, which arc transported 
with x-direction particle velocity. Normally, particles are randomly distributed around the average velocity. 
From statistical mechanics, the moving particles in x-direction can be most favorably described by the 
Maxwellian- Boltzmann distribution function, 

( 2 . 2 ) 

7 T 

where U is the average velocity and A is the normalization factor of the distribution of random velocity, 
which is related to the temperature of the gas flow, i.e. A = — m/2/cT, where m is molecule mass, k is the 
Boltzmann constant, and T is the temperature. 

The transport of any flow quantities is basically due to the movement of particles. With the above 
equilibrium state g , we can split the particles into two groups. One group is moving to the right with u > 0, 
and another group moving to the left with u < 0 . Before splitting the fluxes, let’s first define the useful 
moments of the particle distribution function, 

(' «")„-& = f u n {-) l/2 e^^ du, 

J a ^ 

where integration limit (a, b) of the particle velocity can be (—00, +00), (—00, 0 ) or ( 0 , +00). There is a 
recursive relation for the moments (u n ), which is 

<u" +2 ) a _ 6 = U{u n+ \. b + 

For example, we have 

<«°>„>o = ^erfc( — \/A[/) ; (u °) u<0 = lerfc(VAtf), 

where erfc is the error function, and 

1 p ~ xu1 1 e _At/2 

(u 1 ) u> o = U(u °) u>0 + -^- ; (u 1 ) u<0 = U{u°) u< 0 ---f=-. 

Obviously, if the integration limit is (—00,00), the following relations hold 

<«°> = 1 , (u 1 ) = u. 
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Note that we have dropped the subscript if (a, b) = (—00, 00). 

Depending on the particle moving directions, the total density p can be split into 


pOO 

P + = 9 du 

Jo 

= P{u°)u> 0 


and 


P 



gdu 


= P(«°)u<0- 


Any macroscopic quantities Z without containing explicitly the x-component velocity U , such as density (>. 
y- and z-direction momentum pV and pW , and magnetic field B x B y , can be split si mila rly 

Z+ = Z(u°) u> 0 


and 


Z- = Z{u°) u<0 . 

The above relations mean that the quantity Z is simply advccti :d with the x-direction particle velocity. 
The x— direction momentum pU can be split into 

poo 

(pU)+ — I ugdu 
Jo 

= p{u l ) „>0 

and 

(pU)~ = / ugdu 

J —OO 

= P(U 1 ) U < 0 - 

Similarly, any quantities containing U term, such as B X U, B y U, oU, pVU, pWU , can be split as 

(ZU) + = Z(u l ) u> o 

and 


(zuy = z (u 1 ) u< q. 


For magnetic field, the above splitting implies that the field is frozen into the particle motion and transported 
with the fluid. Note ZU does not include pU , and the splitting of pU will be derived later. 

The splitting of energy can be written as 



= \(n 2 )u>o 
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where pe is the internal energy of the specific distribution function g in Eq.(2.2) with the value of p/4 A. 
Similarly, we have 



= ^ 2 >u<0 

= \pU(u 1 ) u<0 + £^(w°)„< o 
= <« 1 )«<o + pe<u°) u<0 . 


The above equations imply that the kinetic energy \pU 2 can be split as 


-pU 2 


-{\pU 2 ) + + {\pU 2 )~ 

= 2 P U ( ul )u>o + 2 P U ( wl )«<o> 


and the internal energy 


pe = (pe) + + (pe) 

= pe{u°) u> o + pe( u°)„<o- 


In general, besides the thermal energy, pe can include other kinds of internal energy, such as magnetic energy 
in MHD equations. For nonideal gases, the internal energy could have a complicated form as a function of 
p and T. However, the above formulation can still be used to split them in terms of (ii°) u >o and (u°) u< o- 
Since the pressure p is related to the internal energy, it can be split as 

P = P(v°)u> 0 + P(u°)n< 0- 


Now let’s consider the energy transport - J u 3 gdu. The energy transport in the positive ^-direction is 


r°° i 
Jo 2 


7 i u 3 gdu = -{u 3 )„>0 


= {\p U2 + P e )( ul )n> 0 + if/p(u°) u> 0 + ^p(u 1 ) u>0 
= pt(u l )u>0 + lt/p(u°) u > 0 + ip(u 1 )u>0, 

where pe = \pU 2 + pe is the total energy density for the specific distribution g. Similarly, the corresponding 
flux in the negative x-direction is 

rO i i 


/: 


u 3 gdu = -(u 3 ) u < o 


= + pe)(u 1 }„ <0 + ^Up(u°) u<0 + ip(u 1 ) u <o 

= pe(w 1 >u<o + ^Up{u°) u< o + ip{u 1 ) u <o- 


The total energy flux in x- direction is 


/ 


°o i 

1 3 1 tt 2 


2 * u gdu = (- pU 1 + pe)U + pU 
— pell + pU. 
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From the above three equations, we conclude that the total energy flux peU can be split as 

peU = (pcU)+ + {pdjy 

= pe{u 1 ) u >o + pe(u') u< o, 

which is composed of kinetic energy flux splitting 

ipC / 3 = (ipt/ 3 ) + + (ipf/ 3 )- 

= 2 P ^ 2 ( u 1 ) u>o + + 2 P ^ 2 ( u1 )« <0 

and the internal energy flux splitting 


peU = pe(v}) u>0 + pe{u) u<0 . 

At the same time, the splitting of the work done by the pressuie pU term is 
pU = ( P U)+ + ( pU )- 

= ^(f/p(u°) u>0 +p{u 1 ) u>0 ) + ^{Up{u 0 ) u< o + p(u 1 ) u<0 ). 

Note the above splitting formula can be generalized to a hyperbolic system with complicated total energy 
density. 

As a special application of the above splitting principle, let’s split the ID Euler fluxes. The flux function 
for ID Euler equations can be separated into 


/ pU 

pU 2 +p 
\ peU +pU 

where / means free transport. The positive flux Ft is 


= F/fiy 


and the negative part F f is 



P \ 


0 \ 

o 

A 

II 

pU 

+ 

P'W°)«> 0 


j 

P< ) 


|p(“‘)«x + \pU{u°) u>0 ) 

; Ff is 





( P ) 


( ° \ 

Ff = (« x )u<o 

1 pU 

+ 

P<w°>«<0 


\ p< y 


\ 2 p( ul )«<c ^pf/{u°)„< o J 


With the above splitting formula, the numerical flux across a coll interface j + 1/2 for the Euler equations 
can be written as 


F tm = n, + f Pk- 

This is exactly the Kinetic Flux Vector Splitting Scheme for the Euler equations [17, 13], and the positivity 
and entropy condition for the above scheme have been analyzec by many authors, such as [14, 23, 12] and 
references therein. 
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As analyzed in [25], all FVS schemes based on positive (negative) particle velocities suffer from the 
same weakness. The particle free transport across cell interfaces unavoidably introduces large numerical 
dissipation, and the viscosity and heat conduction coefficients are proportional to the CFL time step. In 
order to reduce the over-diffusivity in FVS Schemes, particle collisions have to be added in the transport 
process. 

As a simple particle collisional model, we can imagine that the particles from the left- and right-hand 
sides of a cell interface collapse totally to form an equilibrium state. In order to define the equilibrium state 
at the cell interface, we need first to figure out the corresponding macroscopic quantities qj+ 1/2 there, which 
are the combination of the total mass, momentum and energy of the left and right moving beams. For 
example, for the Euler equations, we have 


Qj+ 1/2 — 



P U 

V PC ' ; + i 

p{u°)u>o 

p(u l )u> 0 

V (pe - lpU 2 )(u°) u> o + \pU(u l ) u> o 
/ p(u°)u<0 

P+ )u<0 

V (pe - \pU’ 2 ){u Q )u< 0 + \pU(v}) u< 0 


j+1 


where pt — \pU 2 is the internal energy density pe. Then, from the “averaged” macroscopic flow quantities 
in the above equation, we can construct the equilibrium flux function 


j+l/2 


pU 

pU 2 + p 
( pc+p)U 


j+ 1/2 


The final flux with the inclusion of both free transport (nonequilibrium) and collision (equilibrium) terms is 


F j+ 1/2 - VFJ+i/2 + (! “ V)Fj+l/2> 

where 77 is a justifiable parameter, which will be analyzed in the next section. The scheme with fixed t] E [ 0 , 1 ] 
is called Partial Thermalized Transport method, which is exactly the first order BGK scheme [25]. With 
the inclusion of equilibrium flux function, the dissipation in the KFVS scheme is reduced substantially. In 
the next section, we are going to extend the above method to MHD equations. In contrast to the Roe’s 
approximate Riemann solver for the Euler equations [18], in the above BGK method, we have strived for even 
less information necessary to form a flux function. So, the above scheme is very efficient. The construction 
of Qj+ 1/2 term at the cell interface gives some ideas about how to construct U± and Mi in the AUSM and 
CUSP-type schemes. It will be interesting to see the results if Ui and Mi are replaced by the equivalent 
values from < 77 + 1/2 term. The splitting of advection and pressure terms in FJ and Fj has the similarity 
with the AUSM-type methods [ 11 , 8 , 20 ]. 
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2.2. Flux Splitting Method for MHD Equations. For the MHD equations, we can use the same 
technique in the last section to split the flux directly. Since the splitting of fluxes is closely related to the 
definition of (u°) and (u 1 ) terms, which are functions of x-diroction velocity U and the “temperature” A. 
For the MHD equations, both gas and magnetic field contribitc to the total pressure p*, and the total 
internal energy is a combination of gas and magnetic energy. With the definition of normal pressure from 
the distribution function g 


/ oo 

(u - Ufgdu = - P -, 

the total pressure (gas + magnetic) in the MHD equations unic uely determines the value of A 

\=J~ = t 

2p* 2p+(fl2 + B 2+ B 2)’ 

where p is the gas pressure. The velocity U in g can be the same as the macroscopic fluid velocity in the 
x- direction. 

After determining A and £/, we are ready to split the MHD flux function, 


F = 


l pU 

pU 2 + Po 

pUV - B x B y 
pUW - B X B Z 
B y U - B X V 
B Z U - B X W 
\ peU + Po U - B x (B y V + B Z W 


= F f + F f 


where Po = p* — B 2 . The positive flux /'V 


f; = (u 1 ) u 


/ 


>0 


+ 


/ p 

pU 

pV 

pW 

By 

B z 

pt 


P o(u°K 


u>0 

—B x B y (u°) u> o 

— B x B z (u°) ti >o 
-B x V(u°) u> o 
-B x W(u°) u> o 
V $(PoU{u°) u >o +Po{u 1 )u>o) - B x (B y V + B z W)(u°) u>0 
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Similarly, the negative flux is 


Ff = (u 1 ),. 


<o 


( P 

P U 
P V 
p\V 

By 

B z 

\ P e 


+ 


P o(u°)u<0 
~B x B y (u°) u< o 
~B x B z (u°) u< o 
-B x V(u°) u< o 

-B x W(u°) u < o 

V ±(p 0 U(u°) u<0 + Po(u 1 )u<o) - + B z W)(u °) u< 0 ) 


Combining the above splitting fluxes, the free transport flux for MHD equations at a cell interface becomes 


pf — p+ _l p~ 

j+ 1/2 “ r jJ ^ r j+i,f 


This formulation is exactly the one given by Croisille et al. [4]. Numerically, the above flux function is very 
reliable and robust [10], and the scheme performs well for these problems for which the Roe scheme fails, 
such as the odd-even decoupling and carbuncle phenomena. However, the accuracy of the above scheme is 
noticeably worse, especially around contact and tangential discontinuities in the MHD applications. 

Now let’s construct the corresponding equilibrium flux for the MHD equations. The corresponding 
macroscopic variables of a equilibrium state at a cell interface arc, 


(2.3) 


where 


Qj+ 1/2 


( p 

pU 

pV 

pW 

By 

B z 

V P* 


\ 


= Qj + Qj+i-i 


) 


3 + 1/2 



< P{u°)v> 0 

pW)u>0 
pV(u°) u> o 

pW(u°) u> o 
By( U °) u>0 

B z (u°) u> o 

V (P e - \pU 2 ){u 0 ) u>0 + \pU{u l ) u>0 


\ 


) 


3 
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and 


1 P{u°)u< 0 

P{U 1 ) u<0 

pV(u°) u< o 

Qj+ 1 = pW(u°) u<0 

By(u°) u<0 
B z ( w °)u<0 

V (pc - \pU 2 )(u °) „<0 + \pb (« 1 > u <o 

With the above “averaged” macroscopic variables q j+ 1 / 2 , the equilibrium flux can be constructed as 

*7+1/2 = * , (?J+l/2 ) 

pU 

pU 2 +p,-B 2 x 
pUV - B x B y 
pUW - B X B Z 
B y U - B X V 
B Z U - B X W 

(pi + p*)U - B X (B X U + B y V + B Z W) 
where B x = B x is a constant in the ID case and 

p. = (7 - 1) (*- - + V 2 + w 2 ) -\(Bl + 1' 2 + B 2 )j + 1 -{B 2 x + B 2 + B 2 ). 

The final flux across a cell interface is a combination of nonequ librium and equilibrium ones, 

(2-4) Fj+ 1/2 = ^^/+l/2 + (1 ~ 7 /) J ^j + l/2» 

where t) is an adaptive parameter. The program from the left and right states to the final flux function is 
given in the Appendix. By removing the contribution from the magnetic field, the above MHD flux function 
reduces exactly to the BGK flux constructed for the Euler equa tions in the last section. 

In the current study, we are more interested in the specific nu nerical flux function for the MHD equations. 
For the lst-order scheme, r/ can be fixed, such as 0.7 or 0.5, in the numerical calculations. Theoretically, the 
parameter r\ should depend on the real flow situations: in the equilibrium and smooth flow regions, the use 
of 7] rsj 0 is physically reasonable, and in discontinuity region, r) should be close to 1 in order to have enough 
numerical dissipation to recover the smooth shock transition. A possible choice for 77 in high-order scheme 
is to design a pressure- based stencil, such as the switch function in the JST scheme [7]. In the high-order 
BGK scheme for the Euler and Navier-Stokes equations [25], wit'i the BGK model as the governing equation, 
the time dependent flux in the gas evolution stage can be obtai ied by following the BGK solution, and the 
relation between the collision time r and viscosity coefficient if well established. For the MHD equations, 
basically we only split the macroscopic flux function without knowing the explicit microscopic transport 
equation for the fluid and magnetic field. Nevertheless, we can ollow the MUSCL-type approach to extend 
the current scheme to high order. For example, we can get the left and right states at a cell interface through 
the nonlinear reconstruction of the initial data, then evaluate the flux according to the formulation given by 
Eq.(2.4). A high-order Runge-Kutta time-stepping scheme is als 3 recommended. For the high-order scheme, 
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the interpolated pressure jump pi and p r around a cell interface can naturally be used as a switch function 
for the parameter ry, such as 


where a can be some constants. 


77 = 1 — exp (—a 


1 Pi - Pr 
Pi +Pr 


), 


3. A Numerical Experiment. For any upwinding schemes, the construction of the flux function, 
or the lst-order scheme, is very important in the understanding of the scheme. Since in the high-order 
extensions, many factors, such as nonlinear limiter, the reconstruction of conservative or primitive variables 
and time* stepping methods, can all affect the performance of the scheme. In the following, we are going to 
apply the current method to the Brio-Wu ID MHD test case [3]. Only the results from first-order method 
with fixed 77 = 0.5 will be presented. 

The initial condition of the Brio-Wu case is 


pi = 1.0 ,£// = 0,pi = l,B Xyi = 0.75, B y j = 1 

on the left, and 

p r = 0.125, Ui = 0,p r = 0.1, B xr = 0.75, B y ^ r — — 1 

on the right. The gas constant 7 is equal to 2 , which corresponds to an internal degree of freedom K = — 1 
for the simulated molecule [25]. Note the gas-kinetic flux splitting formula presented in the last section can 
be applied to any reasonable 7 . 

In order to evaluate the performance of the current method, we are going to compare its numerical 
results with that from the Roe- type MHD Riemann solver [3, 16]. The Roe- type MHD solver is considered 
the most accurate MHD solver existing so far [10], although the robustness of the scheme is questionable in 
some special applications. 

There are 400 grid points used from [— 1 , 1 ] in the ar-direction. The time step is based on At/ Ax = 0 . 2 , 
which is equivalent to CFL number 0.8 in this case. The results at 200 time steps are displayed in Fig. 
4. 1-4.5. The results from the Roe scheme [3, 16], with identical initial condition and time step, arc also 
plotted in these figures. In most regions, the kinetic and Roe- type MHD solvers give almost identical results, 
except the non-conservative quantities at the fast right- moving rarefaction wave. 

Due to the nonconvexity, the MHD equations could present compound waves, which directly connect 
shock and rarefaction. In Table 1 , we list the data at the peak point of the compound wave in the Brio-Wu 
test case. Both results are compared with the theoretical prediction in [3]. Fig. 4.6 gives a close look at 
the density distributions around the right moving shock and the middle contact discontinuity wave. Three 
schemes used here are the current one with 77 = 0.5, Croisille et al.’s KFVS MHD solver, and the Roc type 
MHD solver. The diffusivity of KFVS MHD solver can be clearly observed. 

Table 1 . The Flow Variables at the Peak Point of Compound Wave 


Scheme 

P 

U -Velocity 

V - Velocity 

By 

gas pressure p 

theory [3] 

0.7935 

0.4983 

-1.290 

-0.3073 

0.6687 

Kinetic 

0.8179 

0.4679 

-1.083 

-0.1239 

0.7300 

Roe 

0.8257 

0.4623 

-0.928 

0.0163 

0.7400 
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4. Discussion and Conclusion. In this paper, based on the gas-kinetic theory, we have constructed 
the kinetic flux splitting formula for the MHD equations. We feel that perhaps there is a wide application 
of the splitting techniques presented in this paper. Also, the kinetic flux splitting formulation has the 
similarities with the AUSM and CUSP type schemes [11, 8], vdiere the advection and pressure terms are 
split differently. The numerical results validate the accuracy of the current approach. 

In terms of the current gas-kinetic MHD solver, we have th -3 following remarks: 

(1.) Extension of the current method to the multidimensional case is straight-forward using directionally 
splitting techniques. If there is a jump of magnetic field in the normal direction, such as B x in the x-direction 
across a cell interface, the weakly nonconservative form [16], 


dB x rr dB x 
— - + U — - 
dt dx 


= 0 


can be simply split by changing U in the above equation to U of Eq.(2.3). Also, in order to satisfy V ■ B = 0 
condition, the projection method can be used to clean up the non-zero divergence of the magnetic field [2]. 
(2.) The current scheme is very efficient in comparison with the Roe-type Riemann MHD solver. For 
example, for 1-D calculations, the flux evaluation takes about 1/3 of the CPU time of the Roe-type scheme. 
For 3D calculations, the saving of computational time is enormous. Since we do not use characteristic 
information of the MHD system, the numerical problems related to nonconvexity, non-strictly hyperbolicity, 
and linearization are avoided. Also, the Boltzmann- type scheme is very robust, especially for high-speed low 
density regions [10]. The main reason for this is that the splitting is based on (n n ) u>0 and (u n ) u<0 , which 
accounts for all particle velocities, instead of switching the flix function according to the Mach number 
M > 1 or M < 1 in many other splitting schemes. 

(3.) The extension of the current method to the system with general equation of state p = p(p,e) is 
straightforward. The important point is to distinguish the differences between the splitting of internal 
energy flux peU and the work done by the pressure pU . No singularity and ambiguity in characteristic 
decomposition of the MHD equations will be encountered in the gas- kinetic splitting formulation. 

There arc still many open questions related to the current gas-kinetic approach. First, underlying the 
macroscopic flux splitting, we do not know the exact microscopit equilibrium state for the whole flow system 
including gas and magnetic field. Second, different from the BGK scheme for the Euler and Navier-Stokes 
equations [25], there is no direct way to extend the current metlod to solve dissipative (including resistivity 
and dispersive effects) MHD equations due to the lack of microscopic transport equations, although the 
dissipative terms can be regarded as additional source terms to the current ideal MHD equations. Third, 
in the plasma calculation, particle method is usually used. Hew to make the smooth transition from the 
microscopic particle method to the macroscopic MHD Riemann solver through the gas-kinetic scheme is an 
important and interesting research direction. Even with many unknowns, the potential advantage of the 
kinetic approach over Riemann solver in the construction of namerical flux function becomes clear when 
solving more and more complicated hyperbolic systems. 
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Appendix: Evaluation of Kinetic MHD Flux Functi >n. 

left state = (ADE1 ,AXM1 ,AYM1 ,AZM1 ,AEN1 ,ABX1 ,ABY1, ABZ1) 
right state = ( ADE2 , AXM2 , AYM2 , AZM2 , AEN2 , ABX2 , ABY2 , ABZ2) 
c gas constant GAM $\gamma$, PI=3.14 $\pi$ should ba given. 

c left and right side pressure pi, pr 

APP1= (GAM-1) * (AEN1-0 . 5* (AXM1**2+AYM1**2+AZM1**2) /ADE1) 

* +0 . 5* (2 . 0-GAM) * (ABX1**2+ABY1**2+ABZ1**2) 

APP2= (GAM-1) * (AEN2-0 . 5* (AXM2**2+AYM2**2+AZM2**2) /ADE2) 

* +0 . 5* (2 . 0-GAM) * ( ABX2**2+ABY2**2+ABZ2**2) 

c left and right sides $\lambda$, and macroscopic velocities $U,V,W$ 

AE1=0 . 5*ADE1/APP1 
AU1=AXM1/ADE1 
AV1=AYM1/ADE1 
AW1=AZM1/ADE1 

AE2=0 . 5*ADE2/APP2 
AU2=AXM2/ADE2 
AV2=AYM2/ADE2 
AW2=AZM2/ADE2 

c left and right side particle velocity moments <u"0>, <u~l> 

TEU0=0 . 5*DERFC (-AU1*SQRT(AE1) ) 
TEU1=AU1*TEU0+0.50*EXP(-AE1*AU1*AU1)/SQRT(AE1*PI) 

TGU0=0 . 5* (DERFC (AU2*SQRT(AE2) ) ) 
TGU1=AU2*TGU0-0.50*EXP(-AE2*AU2*AU2)/SQRT(AE2*PI) 

c intermediate (equilibrium ) state, and corresponding $\lambda$ and pressure. 
ADE=ADE 1 *TEU0+ADE2*TGU0 
AU= ( ADE1*TEU1+ADE2*TGU 1 ) /ADE 
AV= ( ADE 1 * AV 1 *TEU0+ADE2*AV2*TGU0 ) / ADE 
AW=(ADE1*AW1*TEU0+ADE2*AW2*TGU0)/ADE 
ABY=ABY 1 *TEU0+ABY2*TGU0 
ABZ=ABZ 1 *TEU0+ABZ2*TGU0 

AE=(AEN1-0.5*ADE1*AU1**2)*TEU0+(AEN2-0.5*ADE2*AU2**2)*TGU0 

* +0 . 5*ADEl*AUl*TEUl+0 . 5*ADE2*AU2*TGU1 
TP= (GAM-1) * (AE-0 . 5*ADE* (AU**2+AV**2+AW**2) ) 

* +0 . 5* (2 . 0-GAM) * (ABX1**2+ABY**2+ABZ**2) 

c gas-kinetic flux function, ETA $\eta$ is a justifiable parameter. 

FM=ETA* (TEU 1 * ADE1+TGU1 * ADE2) + ( 1 -ETA) *ADE*AU 

FU=ETA*(TEUl*AXMl+TGUl*AXM2+(APPl-ABXl**2)*TiU0+(APP2-ABX2**2) *TGU0) 

* +( 1-ETA) *(ADE*AU**2+TP-ABX1**2) 

FV=ETA* (TEU 1 * AYM 1+TGU1* AYM2-ABX1 * ABY 1 *TEU0-A3X2*ABY2*TGU0) 

* +( 1-ETA) * (ADE*AU*AV-ABX1*ABY) 

FW=ETA* (TEU 1 * AZM1+TGU 1 * AZM2-ABX1 * ABZ1 *TEU0-A8X2*ABZ2*TGU0) 

* + ( 1-ETA) * ( ADE*AU*AW-ABX1*ABZ) 

FBY=ETA* (TEU 1 * ABY1+TGU1 * ABY2-ABX1 *AV 1 *TEU0-ABX2*AV2*TGU0) 

* + ( 1 -ETA ) * ( ABY * AU- ABX 1 * AV) 
FBZ=ETA*(TEU1*ABZ1+TGU1*ABZ2-ABX1*AW1*TEU0-ABX2*AW2*TGU0) 

* +( 1-ETA) * (ABZ*AU-ABX1*AW) 

FE=ETA* ( TEU 1 * AEN 1 +TGU 1 * AEN2 

* +0.5* (APP1-ABX1**2) *TEUl+0 . 5*AU1* (APP1- \BX1**2) *TEU0 

* -ABX1* (ABY1*AV1+ABZ1*AW1)*TEU0 

* +0.5* (APP2-ABX2**2) +TGU1+0 . 5*AU2* (APP2- VBX2**2) *TGU0 

* -ABX2* (ABY2*AV2+ABZ2*AW2) *TGU0) 

* + (1-ETA) * ( ( AE+TP) *AU-ABX1* (ABX1*AU+ABY* \V+ABZ*AW) ) 
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Fig. 4.3. y-component velocity distributions with 400 grid points, sol%d line: first-order BGK-type scheme , dashdot line: 
first- order Roe-MHD solver 



Fig. 4.4. B y distributions with 400 gnd points, solid line: first-orde BGK-type scheme, dashdot line: first-order Roe- 
MHD solver 
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Fig. 4.5. Gas pressure p distributions with 400 grid points, solid line: first-order BGK-type scheme, dashdot line: first- 
order Roe-MHD solver 


ngr* moving anocK contact aaoonunuity 




Fig. 4.6. Density profiles around right moving shock and middle contact discontinuity using three first-order schemes, 
+: current kinetic method with r) = 0.5, o; Roe-type MHD solver, *: Croisille et al.’s KFVS MHD solver (corresponding to 
77 = 1.0 in the current scheme). 
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